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Recent neutron scattering studies revealed the three dimensional character of the magnetism 
in the iron pnictides and a strong anisotropy between the exchange perpendicular and parallel 
to the spin stripes. We extend studies of the J\-J2-Jc Heisenberg model with 5 = 1 using self- 
consistent spin-wave theory. A discussion of two scenarios for the instability of the columnar phase 
is provided. The relevance of a biquadratic exchange term between in-plane nearest neighbors is 
discussed. We introduce mean-field decouplings for biquadratic terms using the Dyson-Maleev and 
the Schwinger boson representation. Remarkably their respective mean-field theories do not lead to 
the same results, even at zero temperature. They are gauged in the Neel phase in comparison to 
exact diagonalization and series expansion. The J\-J2-Jc model is analyzed under the infiuence of 
the biquadratic exchange Jbq and a detailed description of the staggered magnetization and of the 
magnetic excitations is given. The biquadratic exchange increases the renormalization of the in-plane 
exchange constants which enhances the anisotropy between the exchange parallel and perpendicular 
to the spin stripes. Applying the model to iron pnictides, it is possible to reproduce the spin-wave 
dispersion for CaFe2As2 in the direction perpendicular to the spin stripes and perpendicular to the 
planes. Discrepancies remain in the direction parallel to the spin stripes which can be resolved by 
passing from S* = 1 to S = 2. In addition, results for the dynamical structure factor within the 
self-consistent spin- wave theory are provided. 

PACS numbers: 75.10.Jm, 75.30.Ds, 75.40.Gb, 75.25.-j 



I. INTRODUCTION 
A. General Context 

Frustrated quantum antifcrromagnctism has been a 
very active field of research over the last 15 years which 
influences many related fields as well^. The frustration 
enhances the importance of quantum effects because clas- 
sical order is suppressed. Hence new and unexpected 
phases may occur and govern the physics at low energies. 
Among the models studied intensely is the J1-J2 model 
on a square lattice^"— and its generalization to three di- 
mensions by stacking planesii"— . Theoretically, two key 
issues arc (i) for which parameters classically ordered 
phases occur and (ii) whether there exists a quantum 
disordered phase between the classically ordered phases. 
The two long-range ordering patterns are either the alter- 
nating order with staggered, Neel type sublattice magne- 
tization or a columnar antiferromagnetic ordering where 
the adjacent spins in one spatial direction (direction a in 
one plane of Fig.[T]) are aligned antiparallel while they arc 
aligned parallel in the other spatial direction (direction b 
in Fig. III). 

The J1-J2 Heisenberg model 

if = Ji^s,.Sj + J2 s^-s,, (1) 

for S ~ 1/2 and its ground states are of broad inter- 
est in solid state physics. The ground state of the sim- 
ple Heisenberg model with J2 = on the square lat- 



tice is the Neel order with staggered magnetization re- 
duced by quantum fluctuationsi^"— . For J2 7^ 0, the 
ground state depends on the ratio Jij J\ of the cou- 
plings. On increasing J2/J1 a value is hit where the 
Neel phase becomes unstable towards a quantum dis- 
ordered state. The intermediate phase is stable in the 
range of 0.4 < J2/J1 ^ 0.6 and is dominated by short- 
range singlet (dimer) formatior^. For J2/J1 > 0.6, the 
spins arrange in a columnar pattern. In the classical limit 
S ^ 00, the transition between the Neel and columnar 
order occurs at J2/J1 ~ 0.5, see Ref. [l9l 

Singh et al^ studied the excitation spectra of the 
columnar phase with series expansion and self-consistent 
spin-wave theory. They calculated the spin-wave veloci- 
ties which depend strongly on the coupling ratio J2I Ji- 
Gapless excitations arc only found dX k = (0, 0) and 
k — (1, 0), while the modes at /c = (0, 1) and k = (1, 1) 
are gapped because of the order by disorder effect. We 
stress that the columnar phase is very well described 
by self-consistent spin-wave theory2iS even in two di- 
mensions and for S = 1/2. The stability of the Neel 
phase, however, is overestimated by self-consistent spin- 
wave theorji^i^ so that the intermediate disordered phase 
is missed. This intermediate phase is seen by a direct 
second order perturbative approacbiS in 1/5' but only 
for spatially anisotropic models with Jia 7^ Jit where 
Jix is the nearest-neighbor exchange in ai-direction with 
X e {o,,b}. The intermediate phase is not the issue of 
the present paper and we point out that it is to be ex- 
pected that it is hardly relevant in three dimensions and 
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for 5 > 



.11.13 



Note that we here and henceforth give all wave vec- 
tors in units of 7r/a where a is the corresponding lattice 
constant. This model was applied to the magnetic exci- 
tations of undopcd iron pnictides^'i^"— i^Sr^. 

It was advocated by Chandra et al^ that stripe order 
can occur at finite temperatures in isotropic frustrated 
Heisenberg models because the stripe order breaks a dis- 
crete symmetry, namely rotation of the lattice by 90°, 
which is not protected by the Mermin- Wagner theorerr^. 
Indeed the transition is of Ising type. This result was 
corroborated by classical^S and semi-classicalSi numer- 
ics. Its significance for the structural and the magnetic 
phase transitions in the pnictides was noted in Refs. [28l 
and[^. 

For completeness, we also mention the additional in- 
stability for ferromagnetic couplings at Ji « — 2J2, 
where the system undergoes a phase transition from the 
coUinear ordered state to a ferromagnetic stato^. 

Since the discovery of superconductivity upon dop- 
ing of the iron-based compound LaOFeAs^. the J1-J2 
Heisenberg model has been used to study the magnetism 
in the parent compounds of the iron pnictide9^ i^^i^^i^°i^^ 
although it does not take the remaining itineracy of 
the charges into account. The magnetic long-range or- 
der is indeed a columnar phase where the spins at the 
iron sites show antiferromagnetic order in a-direction 
and a ferromagnetic order in 6-direction (see Fig. [T]). 
Between the layers (c-direction) , the spins are also 
aligned antiparalleli^"— . The columnar order of the 
spins is also supported by the results of band structure 
calculations2ii2^. 




Figure 1; (Color online) Three dimensional coUinear phase 
with antiferromagnetic interlayer exchange Jc and biquadratic 
exchange Jbq between in-plane nearest neighbor sites. 



A related question still under debate is the reduced size 
of the local magnetic moment. Neutron scattering stud- 
ies find a reduced staggered magnetic moment of (0.31- 
0.41)/iB for LaOFeAs^, (0.8-0.9)/iB for the 122 pnictide 
BaFe2As2^, (0.90-0.98)AiB for SrFeaAsj?^ and 0.8/xb for 
CaFe2As2^i2S. In contrast, theoretical band structure 
calculations determined much higher values, e.g., up to 
2.3nB for LaOFeAs^. 

One conceivable explanation is that the reduction of 



the local moment is due to a strong frustration within 
a localized modelS'i^i'i^'SSi^dS. This requires the sys- 
tem to be in the immediate vicinity of a quantum phase 
transition. But it must be kept in mind that the local 
magnetic moment depends on matrix elements which are 
well known only in the limit of localized electrons. Thus, 
it is possible that the significant reduction of the stag- 
gered magnetization is due to electronic effects such as 
hybridization, spin-orbit coupling and the itineracy of the 
charges^"—. The bottomline of this argument is that the 
value of the staggered magnetization is not a stringent 
constraint for the applicable model for iron pnictides. 



B. Spatial Anisotropy of Exchange Couplings 

Recently, it was shown that a frustrated Heisenberg 
model in the three dimensional columnar phasei^ can 
explain the spin-wave dispersion of CaFe2As2^i^i^ in 
the direction perpendicular to the spin stripes and be- 
tween the layers. However, discrepancies at high ener- 
gies are present in the direction parallel to the stripes. 
These discrepancies can be fitted by a Heisenberg model 
which assumes that the nearest neighbor (NN) coupling 
Ji depends on the spatial direction of the two coupled 
spins, i.e., one introduces Jia » Jn^^nMd^-^. This is 
very remarkable since the difference of the orthorhombic 
distortion between the lattice constants a and b of the 
columnar ly ordered layers is less than l %32,36 ^ which by 
far docs not give reasons for the large spatial anisotropy 
of the magnetic exchange. This argument is further sup- 
ported by the observation that density functional cal- 
culations achieve a good description of the pnictides if 
magnetic columnar order is accounted for. But the con- 
sideration of the orthorhombic distortions is only a minor 
point^. 

Another possible explanation of the spatial anisotropy 
is the possibility of orbital ordering^. But orbital order- 
ing is usually related to higher energies and should lead 
to clear experimental signatures or theoretical signatures 
in density-functional calculations which are so far not 
found. 

So it appears that the magnetism itself generates the 
spatial anisotropy although the original magnetic model 
is not anisotropic. Indeed, the anisotropic order gener- 
ates some anisotropic velocities in self-consistent spin- 
wave theory^iii^. But this effect is not sufficienti^ to ac- 
count for Jia « 40 mcV and Jib w O meV— i^. 

In order to identify a magnetic process which is able 
to generate the observed spatial anisotropy one has to go 
beyond bilinear exchange. This is possible in view of the 
larger spin value S > 1. Indeed, significant higher spin 
exchange processes are inevitable^i^. To be specific, we 
will consider the biquadratic exchange 
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(2) 
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with Jbq > 0. How does a term such as -ffbq influence 
the magnetic excitations? To obtain a rule of thumb we 
adopt a simple view and approximate 

'^bq ^ ^ (Si • Sj) ~ — 2Jbq ^ ^ Si • Sj(Si • Sj) 

+ JbqE(S'-S,)'. (3) 

In the above formula we do not list terms of the types 
sl''^S^/\s^''^sf) or {S^^^sl''>){S^''^sf) because both 
would only contribute for -f ~ S G {x,y,z} and then 
their sum over all spin components would merely yield 
trivial constants. In the direction of alternating spin ori- 
entation one has (Si-Sj) < so that the biquadratic term 
effectively strengthens the bilinear antiferromagnetic ex- 
change: Jf^ > Ji. In contrast, in the direction of par- 
allel spin orientation one has (Si • Sj) > so that the 
biquadratic term effectively weakens the bilinear antifer- 
romagnetic exchange: < Ji. Hence a NN biquadratic 
exchange appears to generate a spatial anisotropy purely 
from the magnetic order. One of the two main goals of 
the present paper is to substantiate this argument phe- 
nomenologically. 

If derived from extended standard Hubbard models 
the biquadratic terms are of higher order in the inter- 
site hoppings than the bilinear exchange coupling which 
implies that it is actually small relative to the bilinear 
exchange^. The exception is a situation where the bi- 
linear terms are subject to near cancellation of antiferro- 
magnetic and ferromagnetic contributions. 

In the pnictides, even the undoped systems are con- 
ducting bad metals. It is known that close to the tran- 
sition from localized to conducting systems higher order 
processes start playing a significant role. For instance, 
the cyclic exchange in one-band Hubbard models be- 
comes as large as 20% of the NN exchange, see for in- 
stance Ref. |50| and references therein. 

In addition, evidence for biquadratic exchange has 
already appeared in calculations based on the self- 
consistent local spin-density approximation (LSDA) by 
Yaresko et al^. Results for the dependence of the 
ground state energy of two intercalated Neel ordered sub- 
lattices on the angle Lp between their sublattice mag- 
netizations are displayed in Fig. [21 With bilinear ex- 
changes no dependence is expected at all. The results 
for BaFe2As2 and LaOFeAs are are reproduced in Fig. 
[2j The dependence of E{ip) = A ■ sit? on </? is com- 
pelling evidence for a biquadratic NN exchange because 
a dependence oc sin^ Lp does not occur in the frustrated 
J1-J2 Heisenberg model. Appropriate values for the bi- 
quadratic exchange are determined from the maximum 
E{'p ^ 90°) in Fig. [2]for 5 = 1. The obtained values are 

Jbq = 21.5 mcV for LaOFeAs, (4a) 
Jbq = 10.1 mcV for BaFe2As2, (4b) 

which are indeed sizeable in view of Ji of the order 40 
meV. 
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Figure 2: (Color online) Dependence of the ground state en- 
ergy on the relative angle tp between the two sublattice mag- 
netizations on the two intercalated Neel lattices, after Ref. 



In general, the electronic situation in the iron pnictides 
is very complex^; up to five bands arc important^i^. 
The appropriate electronic description has not yet been 
identified^i^i^ and thus it is presently not possible to 
discuss the relative size of a biquadratic exchange reli- 
able. Hence we adopt here the phenomenological ap- 
proach to take the microscopic arguments as evidence 
for the existence of such a biquadratic exchange. Next, 
the aim will be to derive the size of Jbq from fits to ex- 
periment. 



C. Methods 

Based on our previous wor h^i^^ , we study the Ji- J2- Jc 
model in the three-dimensional phase with columnar spin 
order. First, we discuss the "critical" value of a; = J1/J2 
where the columnar phase ceases to exist and the cor- 
responding two scenarios for this instability. Further- 
more, we extend the Ji-J2-Jc model by the biquadratic 
exchange discussed above. To this end, an appropriate 
mean-field decoupling has to be introduced which is able 
to tackle biquadratic exchange as well. 

To establish a reliable mean-field approach we employ 
the Dyson-Maleev as well as the Schwinger bosons rep- 
resentation. To gauge the resulting decoupling schemes 
both arc applied to the two dimensional Neel phase of a 
NN bilinear Heisenberg model plus biquadratic exchange 
for S = 1. Then, the successful decoupling is applied to 
the three dimensional columnar phase with biquadratic 
exchange for 1 < 5 < 2. General results are discussed 
before we apply the model to CaFe2As2. 
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II. SCENARIOS FOR THE Ji-J2-Jc MODEL 

Here we discuss two fundamental scenarios for the in- 
stability of spatially anisotropic phases of magnetic long- 
range order. We do not consider first-order transitions to 
other phases; in particular we do not aim to discuss the 
existence of possible intermediate disordered phases be- 
tween Neel and columnar phase. Instead we focus on how 
the columnar phase can become unstable towards fluctu- 
ations. We stress that the self-consistent spin-wave the- 
ory is expected to yield reliable results for S > 1 > 1/2 
and dimension d = 3 > 2 since its results are already 
very good^i^ for S = 1/2 and d = 2. 

It is generally agreed that the staggered magnetization 
in the columnar phase of two-dimensional J1-J2 model 
can take any possible value provided an adequate fine- 
tuning of the coupling ratio x = Ji / J2 is performed2i2^. 
Here we point out that this is no longer true in the pres- 
ence of an inter layer coupling Jc- The quantitative as- 
pects, though not the qualitative ones, depend strongly 
on the spin S. 

Passing from two to three dimensions the additional 
coupling between the columnarly ordered planes cuts off 
the logarithmic divergence of the Goldstone modes^^. 
Thus, in three dimensions the staggered magnetization 
may no longer adopt any arbitrary value > even if the 
ratio X = J1/J2 is increased towards 2. Instead, there 
can be a finite minimum value of the staggered magne- 
tization whose value depends on the relative interlayer 
coupling /i = Jcl J\- 

There are two possible outcomes upon a; — > 2. If the 
value of \x is not too large, it is still possible to drive 
the magnetization to zero while all the three different 
spin-wave velocities remain finite. But for larger values 
of the magnetization remains finite while the small- 
est of the spin- wave velocities vanishesi^. Thus, we face 
two qualitatively different scenarios. They can easily be 
distinguished in plotting the ratio v^/va versus the stag- 
gered magnetization vist- The results for S* = 1/2 and 
S = 1 are shown in Fig. |21 

In the first scenario (solid lines) , the ratio Vb jva always 
stays positive and the instability of the phase is marked 
by the staggered magnetization mst going to zero. Thus, 
we call this scenario the magnetization-driven scenario 
because the breakdown of the columnar phase is driven 
by the vanishing magnetic long-range order. The numeri- 
cal evaluation of the three dimensional integrals becomes 
very difficult close to rrigt = so that the solid curves 
in Fig. [3] can only be computed down to some small fi- 
nite value of the staggered magnetization. Beyond this 
point they are extrapolated as depicted by the dotted 
curves. In this way the phase diagrams of the J\-J2-Jc 
Heisenberg model in Ref. [l^ were determined. 

The dashed lines in Fig.[3]correspond to the second sce- 
nario. We call it the velocity-driven scenario because the 
vanishing of the spin-wave velocity Vh marks the instabil- 
ity of the columnar phase. If x were taken infinitesimally 
higher, would become negative so that no physically 
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Figure 3: (Color online) Ratio Vf^/va as a function of mat 
for S = 1/2 (upper panel) and S = 1 (lower panel). The 
curves either end on the ordinate, i.e, for zero staggered mag- 
netization, or on the abscissa, i.e., for zero spin- wave ve- 
locity < Vc < Va- Thus, two qualitatively different in- 
stabilities emerge naturally. The solid lines correspond to 
the magnetization-driven scenario (vanishing magnetization) 
while the dashed lines correspond to the velocity-driven sce- 
nario (vanishing spin- wave velocity). The dotted lines are 
extrapolations in the magnetization-driven scenario. 



meaningful mean-field model would be found. The stag- 
gered magnetization mgt stops at this point at a finite 
positive value. This finite positive value of the staggered 
magnetization may appear surprising because it signi- 
fies that the magnetic long-range order still persists even 
though the columnar phase ceases to exist because its 
excitations become unstable. 

In the case of S" = 1/2, the transition from the 
magnetization- to the velocity-driven scenario occurs for 
/i,nv ~ 0.134. Thus, the magnetization-driven scenario 
applies in the range < < 0.134. The precise value of 
/i where the magnetization-driven scenario switches over 
to the velocity-driven scenario corresponds to the curve 
in the Vb/va-rUst plane which ends at the origin. If the 
interlayer coupling is increased beyond /imv, the velocity- 
driven scenario applies. 
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For S* = 1, the change from one scenario to the other 
occurs roughly at /Umv ~ 0.01, see lower panel of Fig. 
[3l Thus, the magnetization-driven scenario for 5* = 1 
persists in the range of < /.i < 0.01, which is a much 
smaller region compared to the case for S = 1/2. We 
stress, however, that there is no qualitative difference 
between S = 1/2 and S = 1. For even larger spin 
values the magnetization-driven scenario will only ap- 
ply for a quickly decreasing range of interplane couplings 
ji ~ Jc/Ji- This is due to the fact that the instabil- 
ity of the columnar phase shifts for increasing spin very 
quickly to where it occurs classically^^, i.e., for x = 2. 
Since larger and larger spins correspond to a more and 
more classical behavior it is not astounding that one has 
to draw nearer and nearer to the logarithmic divergence 
of the quantum fluctuations in order for them to become 
important. Hence S oo implies an exponential con- 
vergence a; — 7> 2 and /imv 0. 

We emphasize that the full discussion of the instability 
of the columnar phase requires the consideration of first 
order transitions as well. This is beyond the scope of the 
present paper because we do not know the phase into 
which the columnar phase becomes unstable. It may be 
a disordered phase, but in view of the larger spin 5 > 1 
and the dimension d = 3 we expect that an intermediate 
phase exists only for a very small parameter region. 

From the knowledge for S = 1/2 and 5* = 1 in 
two and three dimensions^"— i^'^i^^i^'^ we presume that 
the magnetization-driven scenarios eventually becomes a 
weak first order transition to an intermediate phase ex- 
isting only within a small parameter region. Further, we 
expect that the velocity-driven scenario becomes a strong 
first order transition to the Neel phase. 



III. BIQUADRATIC EXCHANGE IN THE 2D 
NEEL PHASE 

In this section, we study the effects of a biquadratic ex- 
change (21) on the simple NN Heisenberg model ( J2 = 0) 
with 5 = 1 on the square lattice. The aim is to establish 
a reliable mean-field description. 

We apply the Schwinger bosons as well as the Dyson- 
Malecv representation to the model and introduce the 
corresponding mean-field approximation. Our aim is to 
study the infiuence of the biquadratic exchange on the 
dispersion of the spin-waves. The findings arc checked 
against results obtained by exact diagonalization and se- 
ries expansion. 

The Hamiltonian under study reads 

i? = J ^ S, • S, - Jbq • Sj )' ' 

where J, Jbq > 0. The brackets indicate the sum- 
mation over NN sites. 



A. Dyson-Maleev representation 

First, we apply the Dyson-Maleev transformation^"— 
to the Hamiltonian ([S])- On sublattice A, the expression 
of the spin operators in terms of bosonic operators reads 

S+^bl(2S~h,) (6a) 

Sr = b.^ (6b) 

S^^-S + h^. (6c) 

After a TT-rotation of the spins around the x-axis, the 
transformation on sublattice B is given by 

S+ = hi (7a) 

Sr = (2S-h,)b^ (7b) 

= + (7c) 

The complete Hamiltonian in the Dyson-Malccv rep- 
resentation is given in the Appendix. In preparation of 
the decoupling we introduce the expectation values 

n:={blb^)^(b]b^) (8a) 

a := (bib]) ^ (b^b^) , (8b) 

where i,j are NN sites with i ^ A and j ^ B. All other 
expectation values vanish because of the conservation of 
the total component ^^Sf. For simplicity, we as- 
sume that the expectation values are real. The high-order 
terms are decoupled according to Wick's theorem^. Ne- 
glecting all constant terms the mean-field Hamiltonian is 
given by 

H^'^ = JDM(a) (5 - a) ^ (n, + n, + bib] + b^b^ , 

(9) 

where the subscript DM stands for Dyson-Malccv and 
should not be confused with Dzyaloshinskii-Moriya. We 
define 

JDM(a) := J + [253 - 25^(1 + 5a) 

o - a (10) 

+S'(18a^ + 8a + 1) - 12a^ - 9q;^ - 2a] 

where the new parameter 

a :~ n + a (11) 

has been introduced for convenience. The parameter a 
has to be determined self-consistently, see below. The 
mean-field Hamiltonian ^ can easily be transformed 
into momentum space where it can be diagonalized using 
a standard Bogoliubov transformation. The diagonalized 
Hamiltonian reads 

H^"- E + E^'' ■ (12) 

kGBZ 
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We stress that the full Brillouin zone (BZ), i.e., —it < The constraint 
kry < t: for each component 7 G {a, 6}, is used here and 



hereafter. This is done for simplicity because full trans- 
lational invariance holds and there is only one mode per 
k point; but it does not have a definite value of its total 
Sz component. The dispersion is given by 



a]a, + b]b, = 2S 



Wk = 2JDM(a) (S - a) 
Wk = 4 — (cos ka + COS kbf' 

and the ground state energy by 



E^^ = JDM(a) {S-a)E 



MF 



E 



MF 



E ('^k - 2) , 



(13a) 
(13b) 



(14a) 
(14b) 



keBZ 



The gapless excitations at fc = (1, 1) which is the mag- 
netic ordering vector of the Neel phase and at k = (0,0) 
are the expected Goldstone modes, because the ground 
state of the system has broken symmetry^. The exis- 
tence of these gapless excitations is guaranteed by the 
systematic expansion in 1/S. We draw the reader's at- 
tention to the fact that this argument ensures massless 
modes only for the systematic expansion, i.e., for a par- 
ticular value of a. But it turns out that the vanishing 
of the energy of the Goldstone modes does not depend 
on the precise numbers of the expectation values so that 
it does not matter whether the expansion is performed 
systematically or self-consistently. 

The self-consistent equation for the parameter a is 
found by comparing ^ with the mean-field ground state 
energy (fTH) per site yielding a = (l/2)£;^^/7V if N is the 
number of sites. Hence we have in the thermodynamic 
limit N ^ 00 



1 1 



a 



d^k (wk - 2) 



BZ 

= -0.0789737105. 



(15a) 
(15b) 



One integration in (jl5bp can be done analytically, the 
second one with any computer algebra system. Note that 
due to the simplicity of the system no real self-consistency 
needs to be determined; a can be computed directly. 



B. Schwinger boson representation 



(17) 



restricts the infinite bosonic Hilbert space to the rele- 
vant physical subspace of the spin S. In this way, the 
Hamiltonian is rewritten as 



4^y- 



(18) 



where the bond operators 



4j = ala] + bib] 



A, 



aid J +Wbj 



(19a) 
(19b) 



are used. They connect adjacent sites on the different 
sublattices with i E A and j <E B. 

The mean-field approximation in terms of the bond 
operators is based on an expansion in the inverse num- 
ber of boson flavors 1 /A/^i^^i^. Thus, we intermediately 
extend to J\f Schwinger boson flavors to justify the ap- 
proximation. For J\f flavors, the bilinear term reads 



with the generalized bond operators 



A] 



AT 

E °l 

m— 1 



A. 



E« 

m— 1 



and the constraint 



E 



m— 1 



We define the expectation value 

A■.^{A^^^{A\^ >o. 



(20) 

(21a) 
(21b) 

(22) 
(23) 



which is proportional to M according to the generalized 
constraint (22). Hence, the bilinear term is decoupled 
according to 



Here we apply the Schwinger boson representatior^i^ to 
the Hamiltonian ([5]). The spin operators are expressed 
as 



St 

s- 



(16a) 
(16b) 

(16c) 



JiMF 



h. c. 



(24) 



where we only keep the leading order O ((l/A/")*^) omit- 
ting constants. 

In complete analogy, the quartic term of bond opera- 
tors is decoupled in leading order by 



jp (4^»j) ■ (44j) ~ -jp {a\^a + aa^^ 



(25) 
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The factor 2 is a consequence of Wicks's theorem, as there 
are two possibihties to contract the operators or A^j . 
Thus, the mean-field decoupling of the entire biquadratic 
term reads 



MF 



2^ 



A 1 



A^ 



AfS- 



/4t 



A- 



(26) 



Since A oc AfS, A? > MS"^ holds so that the biquadratic 
term yields a positive contribution in total. Due to the 
prcfactor — Jbq, cf. Eq. (1181) . the biquadratic term en- 
hances the bilinear one 
Returning to A/" = 2 
finally given by 



the mean-field Hamiltonian is 



MF 



H 



where 



A] 



(27) 



B := AJsb{A)/2 (28) 
Jsb{A) :=J- 2Jbq52 (1 _ Ay{2S^)) . (29) 

The Lagrange multiplier A is introduced to enforce the 
constraint (jl7[) on average. In the symmetry broken 
phase at zero temperature, A is fixed to the value 



X = AB = 2 JsB (A) A 



(30) 



in order to retrieve massless Goldstone bosons. In mo- 
mentum space, the mean-field Hamiltonian (j27p is diago- 
nalizcd by a standard Bogoliubov transformation leading 
to the dispersion and the ground state energy 



^^k = Y -^^ ^ (4S (cos ka + cos fcf,)) 



(31a) 
(31b) 
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The self-consistency equations for the parameters A 
and A arc given by 



25: 

-4A 



1 dE^^ 
N dX 
1 dE^^ 
N dB 



(32a) 
(32b) 



These two equations can be combined to yield A in the 
symmetry broken phase 



2ii (2^)' 

BZ 



\j ^— (cos ka + COS kbf (33) 



in the thermodynamic limit. Note that the massless 
Goldstone modes are again guaranteed by the system- 
atic expansion which is here done in the inverse number 
of flavors 1 /TV. As for the Dyson-Maleev mean-field ap- 
proach the expansion can also be done self-consistently 



without spoiling the Goldstone theorem because the pre- 
cise number of the expectation value A does not matter 
for this qualitative aspect. For S* = 1, the calculated 
value is 



A = 2.1579474210 



(34) 



which corresponds to 2(1 — a) from the Dyson-Maleev 
calculation. Thus the obtained dispersions are identical 
for zero biquadratic exchange because then JgB = Jdm 
holds. 



C. Results 

We are interested in the influence of the biquadratic 
exchange on the excitation energies. Thus the ratios of 
the dispersions with and without biquadratic exchange is 
an appropriate quantity. Since the shape of the disper- 
sions is the same for the Dyson-Maleev and the Schwingcr 
boson representation, no more quantities need to be com- 
pared. This ratio depends linearly on the coupling ratio 
V Jhq/J, see Eqs. pU]) and (^5]) . This behavior is de- 
picted in Fig. m and the corresponding slopes are given 
in Tab. D 




V=Jbq/J 



Figure 4: (Color online) Influence of the biquadratic exchange 
on the excitation energies in the 2D Neel phase {S = 1). The 
dashed lines are fits to the corresponding data points. Note 
the almost perfect straight lines obtained in exact diagonal- 
ization and in series expansion. 

The mean-field results for the 2D Neel phase with bi- 
quadratic exchange arc checked to the results obtained by 
other methods. In detail, the excitation energies can be 
calculated by an exact diagonalization of the Hamiltonian 
([5l) on a finite 4x4 lattice^ and by scries expansio n^^-^'^ . 
For both methods, the slope is determined by a linear 
fit, see Fig. H] All results are compared in Tab. U and 
displayed in Fig. SI The data points are the maximum 
excitation energies normalized to the maximum without 
biquadratic exchange. 
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lilt! LllUU. 




Schwinger bosons 


2.65674 


Dyson-Maleev 


1.27708 


exact diagonalization (4x4) 


1.24505 


series expansion 


1.13423 



Table I: (Color online) Slopes for the influence of the bi- 
quadratic exchange on the maximum excitation energy rel- 
ative to the one without biquadratic exchange, see Fig. 3] for 
the data from which the slopes are fitted for exact diagonal- 
ization and series expansion. 



The Schwinger bosons result immediately catches one's 
eye because the slope is far too large by a factor greater 
than two. In contrast, the Dyson-Maleev result matches 
very well with the exact diagonalization. Of course, the 
exact diagonalization data is hampered by finite size ef- 
fects even though they are only moderate because we 
consider the ratio of the maximum of the dispersion. The 
slope of the series expansion, which is essentially exact, 
is slightly smaller. The deviation of the Dyson-Maleev 
slope from the series expansion is in the range of 10- 
15%, which is well acceptable for a mean-field approach 
applied to complicated terms made from eight bosonic 
operators. Extending the model to three dimensions, we 
expect the result to improve because mean-field approx- 
imations generally become the more accurate the higher 
the number of dimensions. Therefore, we choose to em- 
ploy to the Dyson-Maleev representation and the corre- 
sponding mean-field approximation in the study of the 
columnar phase. 

Discrepancies of the Schwinger boson mean-field the- 
ory to other treatments appeared already without bi- 
quadratic exchange when the theory was introduced by 
Auerbach and Arovas^^iffi. In their calculations, the 
mean-field free energy exceeded previous results by a fac- 
tor of 2 and the sum rule of the dynamic structure factor 
exceeded its exact value by a factor of 3/2. 

In the mean-field approximation for Schwinger bosons, 
we treat both Schwinger bosons on a site as independent. 
Thus, the constraint p7)) is violated because a change in 
the occupation of the a boson should always be connected 
to a change of the occupation of the b boson on the same 
lattice site. Hence, we suggest that the additional factors 
appear because of the missed suppression of fluctuations 
of the boson number due to the constraint. Since the 
biquadratic exchange can be roughly viewed as a bilinear 
exchange multiplied by a NN expectation value, see Eq. 
([3]), it is comprehensible that the slope of the Schwinger 
boson mean-field result is too large by about a factor of 
2 because the NN expectation value is overestimated by 
this factor—. 



IV. 3D COLUMNAR PHASE WITH 
BIQUADRATIC EXCHANGE 

As explained in the Introduction we advocate that a 
localized spin model for the pnictides has to comprise 
a biquadratic term in order to account for the spatial 
anisotropy of the spin-wave velocities measured by in- 
elastic neutron scattering (INS) and computed by den- 
sity functional theory. Thus, we set out to investigate 
the Hamiltonian 

H = Ji Si ■ Sj + J2 Si • Sj 

+ t/c ^ ^ Si ■ Sj — </bq ^ ^ (Si ■ Sj) , 

[*>i] (iJ) 

where a NN in-plane biquadratic exchange {Jhq > 0) is 
introduced. The brackets and correspond to 

in-plane nearest and next-nearest neighbor sites, while 
the bracket indicates the exchange between inter- 
plane nearest neighbor sites. The very small orthorhom- 
bic distortion in the columnarly ordered layers is ne- 
glected. Justified by the checks in the previous Section, 
we employ the Dyson-Maleev representation and the en- 
suing mean-field approximation. In order to establish our 
notation we provide a brief derivation of the spin-wave 
dispersion and the self-consistency equations. For further 
details the reader is referred to Ref. [6J. 

For the mean- field decoupling, the following parame- 
ters are needed: 

• average occupation number per lattice site 

^ ■■= (bib.) = {b]b^) ^ (36a) 

• in-plane NN antiparallel spin orientation perpen- 
dicular to the spin stripes (a-direction) 

«i := (blbl) = {b,b^) (36b) 

• interlayer NN antiparallel spin orientation (c- 
direction) 

«c {bjb]) - (b^b^) (36c) 

• in-plane next-nearest neighbor (NNN) antiparallel 
spin orientation 

02 (blbl) = {b,b^) (36d) 

• in-plane NN parallel spin orientation parallel to the 
spin stripes (6-direction) 

f ■■= (bh) = (bh) (36e) 
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:= n 


+ fli 


etc 


:— n 


+ flc 


a2 


:= n 


+ 02 


13 


n 


-/• 



It turns out to be advantageous to introduce the com- 
bined parameters ai, 012, etc and /? according to 



(37a) 
(37b) 
(37c) 
(37d) 



The values of these parameters are determined by the 
self-consistency equations given below. 

Thereby, the mean-field decoupling of the Hamiltonian 
?5l) reads 



H = H^+H\\+ i?e + i?NNN (38a) 
Hi_ = J2Xia {S - ai) ('\ + + ^Ib] + b-b 

(38b) 

= -J2XU {S~/3)J2 + ^ - ~ b]b,) 

(38c) 

He = J2Xfi (S - ac) (n^ + + bib] + b^b 

(38d) 

i^NNN = ^2 (S* - 02) Y + + ^1^1 + h bj^ 

(38e) 

with X := J1/J2, n ;= Jc/Ji, v := Jbq/Ji and 

xia x + [2S'3 - 2S'2(1 + 5ai) + 

D — q;i 

^(ISa? -h 8ai + 1) - 12a? - 9aJ - 2ai] (39a) 



Xib := x 



XV 



[2S^ - 252 (1 + 5;3) 



5(18^2 _ ;L2^3 _ g^2 _ ^ 



(39b) 



The mean-field Hamiltonian (p8| is easily transformed 
into momentum space and diagonalized by a standard 
Bogoliubov transformation. Thereby, the spin-wave dis- 
persion 



(40) 



is obtained with 



Ak := J2[2(S' - as) + xia{S - ai) + 

x^i{S - ac) + xib{S - l3){cos h - 1)] (41a) 
-Bk := J2[xia{S — ai) cos ka + x^{S — ac) cos kc+ 

2(5* — 02) cos fca cosfcfe]. (41b) 

The components ka, kf, and kc of the momentum vector 
k are directed along the crystal axes as shown in Fig. [T] 
The dispersion is gapless at fc = (0, 0, 0) and k = (1, 0, 1) 
corresponding to the required Goldstone modes^iii^. Note 



that this feature is again guaranteed by the systematic 
expansion in the inverse spin 1/5. As for the Neel phase 
the expansion can also be done self-consistently without 
spoiling the Goldstone theorem because the precise num- 
bers of the expectation values do not matter for this qual- 
itative aspect. But we consider this a highly non-trivial 
aspect in view of the four different quantum corrections 
occuring here. 

For small momenta, the dispersion (|40)) can be ex- 
panded 



« ^Jvlkl+vlk^+v^ck'c (42) 
and one obtains the spin-wave velocities 

vl = i2J2)^[2{S-a2)+xiaiS-ai)+xfi{S-ac)]x 
[2(8 ~ a2) + xiaiS ~ ai)] (43a) 

vl = {2J2)^[2{S^a2)+xiaiS-ai)+xfiiS-ac)]x 
[2(5-a2)-xib(5-/3)] (43b) 

vl = i2J2)^[2iS - a2) + xiaiS ~ ai) + x^i{S - ac)]x 
xfi{S-ac). (43c) 

The parameters of the quantum corrections are deter- 
mined from the self-consistency equations in the thermo- 
dynamic limit 



1 1 



(2^)^ 



d-^k 



Ah - Bv cos kn 1 



BZ 



1 1 



a2 



2 {2^f 

1 1 
2(^ 



BZ 



A^k 



d-^/c 



Ak — B\f^ cos kc 1 



(44a) 



— By cos ka COS fcfo 1 



(44b) 



BZ 



= - 



1 1 



2 (2^)' 



BZ 



(44c) 



where the integrations run over the full Brillouin zone 
(BZ) as before. The above set of equations is solved 
by numerical integration and by four-dimensional root 
finding. 

The staggered magnetization mgt is defined as nist := 
(Sf) (-1)'^ where cr = for i e A and cr = 1 for i e B. 
In the thermodynamic limit, it reads 



TOst = S" + 



1 1 1 



2 2 (2^)^ 



BZ 



(45) 



A. General Results for S = 1 

We discuss the qualitative aspects for 5 = 1 and 
^ = 0.25, which is a generic value for the relative in- 
terlayer coupling. The values of the relative biquadratic 
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Figure 5: (Color online) Staggered magnetization as a func- 
tion of X for S — 1 and fi = 0.25. All curves are in the 
velocity-driven scenario which we presume to signal a strong 
first order transition to the Neel phase. The magnetization re- 
mains finite at the endpoint where the columnar phase ceases 
to exist, for details see Sect. HIl 



exchange f are chosen in the range of 0.1 to 0.7 as mo- 
tivated in Sect. II Bl For comparison, wc also show the 
results without biquadratic exchange. 

From the staggered magnetization in Fig. [SJ we con- 
clude that the biquadratic exchange destabilizes the 
columnar phase and drives the critical point towards 
lower values of a; = Ji/J2- There is only a negligible 
influence on the maximum value of rrist, which is also 
shifted further left. 

Important is the effect of the biquadratic exchange on 
the renormalization of the quantum correction parame- 
ters (for definitions see Eqs. ([371 EH EH) ) shown in Fig. [HI 
The effective exchange perpendicular to the spin stripes 
is significantly strengthened by factors up to 2, sec upper 
left panel of Fig. [51 In contrast, the effective exchange 
parallel to the spin stripes stays almost constant except 
for small frustrations and in vicinity of the end point, see 
upper right panel of Fig. [5] Because the biquadratic ex- 
change in the Hamiltonian psp is restricted to in-plane 
NN sites, the only effect on the renormalization of the 
interlayer and NNN exchange is the shift of the critical 
point to smaller values of a;, see lower panels of Fig. [HI 
Hence, wc indeed find the expected effect that an NN bi- 
quadratic exchange enhances the spatial anisotropy. This 
validates our qualitative considerations in Sect. II Bl 

The ratios of the spin- wave velocities depicted in Fig. [7] 
are weakened by increasing biquadratic exchange because 
the strengthening of the effective exchange perpendicu- 
lar to the stripes leads to a greater spin- wave velocity Va 
and thus to a smaller ratio Vb/va- In general, the qualita- 
tive behavior of the ratios is similar to the one observed 
without biquadratic exchange. 

All in all, a biquadratic exchange acting on in-plane 
NN sites increases the anisotropy of the exchange par- 
allel and perpendicular to the spin stripes. The grow- 



ing anisotropy with increasing biquadratic exchange v = 
Jhq/ Ji is caused by the strong renormalization of the 
exchange perpendicular to the spin stripes. The ex- 
change parallel to the spin stripes experiences a marginal 
strengthening and is weakened only in proximity to 
the critical point. Hence, even a large biquadratic ex- 
change does not induce an effective ferromagnetic ex- 
change J^^ ^ in 6-direction. 

To render this point quantitatively we define the effec- 
tive exchange couplings 



Teff 

Teff 
•Jib 

reff 



-off 



= J22:ia(l - Cei/S) 
= J2Xlbil-f3/S) 

^ Ml - aJS) 
= J2(l - a2/S) 



(46a) 
(46b) 
(46c) 
(46d) 



according to the mean-field Hamiltonians (j38bl I38cl I38dl 
These definitions enable a direct comparison of 



the effective spatial anisotropy to the ones determined 
in fits based on linear spin-wave theory as they arc used 
in experimen t'^^1^'^ , for the analysis of density-functional 
theory^ and in other theoretical analyse a^^i^^ The rela- 
tive anisotropy is shown in Fig. |S] for the two-dimensional 
model, i.e., for Jc = 0. The results for finite three- 
dimensional coupling are qualitatively the same. 

Strikingly, the size of the spin really matters due to 
the different relative importance of quantum fluctuations. 
For S — 1 the self-consistency procedure prevents the 
occurrence of negative effective couplings in b direction 
- even for very large biquadratic exchange. This can 
be understood by inspecting Eq. (j39bp . For large u the 
expectation value /? goes to zero. Hence the influence 
of the square bracket multiplying v decreases more and 
more because the terms independent of /3, i.e., 25'^ — 25^, 
cancel for 5 = 1. Hence no zero or negative effective 
coupling along the spin stripes occurs. Only for S > 
1 a zero or even negative effective coupling along the 
direction of parallel spins is possible. We will come back 
to this point when comparing with experimental data. 



V. APPLICATION TO THE IRON PNICTIDES 

Here we discuss the applicability of our model to iron 
pnictides. We consider S' = l,5' = 3/2 and S = 2. 
The former value is suggested by the relatively low values 
for the staggered magnetizations measured^iiSi^^i^S, see 
also list in Ref. [l^. Further support for this value comes 
from the successful studies of two-band models^"—. 
Another interesting support is provided by advanced 
Gutzwillcr calculations of the distribution of local charges 
in five-band models which are consistent with S = 1—. 
On the other hand, simple chemical electron count im- 
plies that the iron is doubly positive charged so that the 
d-shcU contains four holes and Hund's rule implies S = 2. 
Indeed, recent findings for closely related iron compounds 
showed that static local moments up to 2.2 — 3.3/iB can 
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Figure 6: (Color online) Renormalization of the effective exchange parameters by quantum fluctuations as function of x for 
5=1 and jj, = 0.25. All curves are in the velocity-driven scenario. 



occur—"— or even larger (4.7/iB)^- Furthermore, we 
stress that the value of the staggered static magnetiza- 
tion per site represents only a lower bound for the local 
moment. Moreover, the observable static moment also 
depends on electronic matrix elements^ influenced by the 
complicated electronic situation. In view of these ambi- 
guities, it is certainly worthwhile to consider also S* = 3/2 
and 5 = 2. 

Experimentally, the magnetic dispersion is best known 
for the 122 compound CaFe2As2^i^i^. Hence we fit the 
dispersion obtained by self-consistent mean-field theory 
of the three-dimensional model p5[) with biquadratic ex- 
change to the measured spin-wave dispersion. Then we 
are able to draw conclusions regarding the quality of the 
agreement and the plausibility of the exchange values ob- 
tained. 

The results of fits to Zhao's data are given in Tab.HIlfor 
plausible values of Jbq. In order to fix Jbq independently 
an additional fourth piece of information is required from 
experiment, for instance the energy at k = (0, 1, 0). This 
energy was only determined by Zhao et al^. The curves 
shown in Fig. [9] illustrate that perfect agreement is pos- 
sible for small energies and perpendicular to the spin 



V 


0.3 


0.4 


0.5 


0.6 


3.0 


X 


0.645 
0.297 


0.616 
0.314 


0.589 
0.332 


0.565 
0.349 


0.284 
0.793 


Jl 

Jc 
■h 
Jbq 


18.9 
5.6 
29.4 
5.7 


17.9 
5.6 
29.0 
7.1 


16.9 
5.6 
28.7 
8.5 


16.0 
5.6 
28.4 
9.6 


7.1 
5.6 
24.9 
21.2 


Tott 

Jla 

TOff 

''lb 
T-efl 

Tcfl 

^2 


25.4 
18.8 
5.3 
31.1 


26.3 
17.9 
5.3 
30.7 


27.1 
17.0 
5.3 
30.3 


27.9 
16.3 
5.3 
29.9 


36.2 
8.0 
5.3 
25.8 



Table II: Fit parameters of the model (|35p in the three di- 
mensional columnar phase for given values of the relative 
biquadratic exchange v = Jh^/Ji- The parameters are de- 
termined by fits to the experimental spin-wave velocities in 
CaFe2As2^. The exchange constants J°^ are the exchange 
constants which would provide the same results for a model 
without biquadratic exchange in linear spin- wave theory, see 
the definitions (|46p . 



stripes, i.e., the direction of parallel spins. We refrain 
from showing the results for Diallo's experimental data. 
The results are very similar, despite the smaller interlayer 
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X=Ji/J2 V=Jtiq/J^ 



0.4 




X=J.|/J2 

Figure 7: (Color online) Ratios of the spin-wave velocities as 
a function of x for S = 1 and = 0.25. All curves are in 
the velocity-driven scenario, i.e., the ratio vt/va vanishes at 
the end point. This is at the limit of the resolution of the 
numerical evaluation of the three-dimensional integrals. 



coupling derived from Diallo's INS data. The interested 
reader is referred to Ref. [6^ 

But even for an artificially large biquadratic exchange 
v = 3 the qualitative behavior of the dispersion along 
the spin stripes is not described satisfactorily, see right 
panel of Fig. [HI At first glance, this comes as a surprise 
since naively a sufficiently large biquadratic exchange 
should induce an arbitrary spatial anisotropy. But the 
self-consistency prevents this, see Fig. \E\ and Eq. (j39b|) 
and the discussion of them. 

The understanding of the quantitative effect of the bi- 
quadratic term tells us that only a larger spin value may 
lead to the experimentally observed anisotropy. Indeed, 
for S = 2 and for S ~ 3/2 we achieve a perfect fit, see 
Tab. mil But for S" = 3/2 the required value of Jbq ap- 
pears unreasonably large relative to Ji. The resulting 
dispersion is shown in Fig. [TU] for 5 = 2. The disper- 
sion obtained for 5 = 3/2 (not shown) look essentially 
the same. Note the excellent agreement obtained with- 
out assuming any spatial anisotropy in the model itself. 
It is the magnetic long-range directional Ising-type order 



Figure 8; (Color online) Relative anisotropy, cf. definitions 
(|46p . as function of the relative biquadratic exchange v in the 
two-dimensional model for S = 1, 3/2 and S — 2. 
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TCfl 
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icfl 
Jib 


TCfl 


icff 


3/2 


0.52 


0.56 


1.20 


12.4 


50.4 


-5.7 


5.2 


18.8 


2 


0.78 


0.37 


0.37 


9.3 


50.4 


-5.7 


5.2 


18.8 


Zhao et al. 










49.9 


-5.7 


5.3 


18.9 



Table III: Fit parameters of the model (|35|l in the three di- 
mensional columnar phase for S = 'i/2 and S = 2. The 
parameters are determined by fits to the experimental disper- 
sions in CaFe2As2^. Exchange couplings in meV. For com- 
parison, the exchange constants Jf^ obtained in linear spin 
wave theory by Zhao et alr^ are given. 

which induces the strong spatial anisotropy. We judge 
the fit parameters necessary for 5 = 2 to be perfectly 
reasonable, in particular the moderate value of the bi- 
quadratic exchange of 37%. 

The above finding provides an interesting piece of in- 
formation for the description of the magnetic excitations 
in the pnictides by a model of localized spins. But it 
leaves open the issue why the observed staggered mo- 
ments are much lower than 4/iB which is the value one 
would expect for 5 = 2. Here we can only speculate that 
the complicated local electronic levels and the remaining 
itinerant character of the charges are the physical reasons 
for this. Also, the issue of line broadening due to Landau 
damping is not included in the present mode P^i''^'^^ . 

All in all, the additional biquadratic exchange influ- 
ences the dispersion parallel to the spin stripes and 
strengthens the anisotropy of the effective in-plane NN 
exchange constants. However, for 5 = 1 it is not possible 
to reproduce the whole experimental dispersion measured 
by Zhao et al. This is possible for 5 = 3/2 for parame- 
ters which do not appear to be realistic. For 5 = 2, the 
dispersion can be described very well for realistic param- 
eters. The question why the observed magnetic moments 
are much smaller than one would expect for 5 = 2 re- 
mains unresolved at present. The local electronic situ- 
ation and the residual itineracy are likely candidates to 
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Figure 9; (Color online) Spin-wave dispersion ()40|l resulting from the self-consistent mean-field theory of the model (|35p for 
5 = 1 in the three dimensional columnar phase with biquadratic exchange for various values oi v = Jhq/Ji- The dispersions 
(lines) are plotted for the parameters given in Tab. [Ill The red dots are experimental data extracted from INS for CaFe2As2 
from Ref. [33 . Wave vectors are given in units of tt over lattice constants assuming an orthorhombic crystal. Note that along 
the path in the BZ shown in the left panel all fits collapse in one curve. Significant differences occur only along the path shown 
in the right panel. 
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k=(ka,0,0) k=(0,0,y k=(0,\<^,0) 

Figure 10: (Color online) Same as in Fig. [9] but for S = 3/2 and 5 = 2, the parameters are given in Tab. IIIll 



reduce the magnetic moment. 

VI. DYNAMIC STRUCTURE FACTOR 

Besides the dispersion it is the dynamic structure fac- 
tor S'(k, w) which matters for the understanding of ex- 
perimental observations. In addition to the information 
about the energies of the collective excitations, S'(k,cj) 
contains information about the relevant matrix elements. 
In the Dyson-Malecv representation, the inelastic part of 
the dynamical structure factor for spins at T = reads 

(k, = NniS- n) \ 5 (u; - u;^) , (47) 

where n, and i?k are defined in Sect. IIVI Because of 
the rotational symmetry about , ^^^(k, co) is identical 



to S'j?{k^ uj). In the limit of k ^ (0, 0, 0), the dynamical 
structure factor (|47)) vanishes because = -Bk- In the 
limit of k — >• (1,0,1), the dynamical structure factor (|T7)) 
diverges because Ak = — Bk- Note that in both cases we 
have Wk — ^ 0. Similar results have been derived within 
linear spin- wave theory in Ref. IT^ 

Constant-energy cuts are computed from (|T7)) for 
equally weighted twinned domains. This means that the 
spatially anisotropic result from (|47p is superposed for 
one choice of the crystallographic a and b directions and 
for the swapped choice a •H- & so that the resulting su- 
perposition is spatially isotropic in the sense that there 
is no difference between the a and the b direction. The 
results for S = 1 are shown in Fig. [TTJ those for 5' = 2 in 
Fig. [T^l They are to be compared with the experimental 
findings in Ref. IH. 

For low energies, concentric rings emerge from the 
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{e)E = 137 ± 15 mcV, kc = 4, (f )£ = 135 ± 10 meV, = 4.5 {g)E = 144 ± 15 meV, fcc = 5 {h)E = 175 ± 15 mcV, fcc = 5.2 

Figure 11; (Color online) Constant-energy cuts of the dynamical structure factor per spin (|47|l integrated over the given 
energy interval for a twinned sample for 5=1 and the relative biquadratic exchange v = 0.5. Other parameters are given in 
Tab. ini The range of ka^b and the values of fee and E have been chosen to match the rendering of the experimental datai^. 
The reciprocal lattice vectors are given in units 7r/lattice constant. Note that the weights in the given energy intervals are 
dimensionless according to formula ((T7)l. 



magnetic ordering vector Q — (1,0,3). They display a 
certain ellipticity which is not surprising in view of the 
spatial anisotropy of the spin order. The rings increase 
in size for higher energies. For 5* = 1 (Fig. [TT|) . less in- 
tensive spots (1,1, fee) appear additionally for E = 100 
meV and = 115 meV which merge with the concen- 
tric spin- wave rings for higher energies. For S ~ 2 (Fig. 
I12p . the circular signatures of the dispersion cones persist 
up to 144meV; only for even higher energies significant 
additional features occur. In general, one observes the 
trend of decreasing intensity for increasing energy be- 
cause 5'(k, w) cx 1/cJk, see Eq. (|47| . 

The comparison to the experimental datai^ shows that 
the S = 2 results match better which is not surprising 
since they reproduce the dispersion everywhere, see Fig. 
\TU[ The experimental data displays also rings of high 
scattering intensities which increase very much like the 
theoretical results do. But the experimental results are 



broader indicating a finite life-time of the magnetic ex- 
citations. This effect is lacking in our model for two 
reasons: (i) The spin-only model does not include any 
Landau damping due to the decay of magnons into elec- 
tronic particle-hole pairs. The consideration of this effect 
would require to pass to a doped t-J type of model or to 
switch completely to an itinerant approach, see for in- 
stance Ref. l72i . (ii) Even in the framework of a spin-only 
model the scattering of magnons from other magnons will 
lead to life-time effects which are not taken into account 
by a mean- field approach. 



Hence, the agreement of Fig. [12] to the experimental 
data is encouraging, in particular in view of the neglect of 
damping effects. We conclude that the description of the 
magnetic excitations of undoped pnictides within a model 
of localized spins is possible if a significant biquadratic 
exchange is included and = 2 is chosen. 




0.5 1 1.5 0.5 1 1.5 




{e)E = 137 ± 15 mcV, fc^ = 4 {i)E = 135 ± 10 meV, fc^ = 4.5 



Figure 12: Same as in the previous Fig. 1111 but 

VII. CONCLUSIONS 

In this paper, we obtained the following results: 
First, we pointed out that the instability of the colum- 
nar phase, i.e., the phase with spin stripes, can be- 
come unstable in two qualitatively different ways in three 
dimensions. Either the sublattice magnetization van- 
ishes so that the long-range order parameter vanishes 
(magnetization-driven scenario). Or the smallest spin- 
wave velocity vanishes (velocity-driven scenario). The 
latter scenario is not possible in two dimensions because 
a vanishing velocity implies a logarithmically diverging 
correction of the quantum corrections of the magnetiza- 
tion so that the magnetization will always vanish before 
the velocity becomes zero. But the subtlety of a logarith- 
mic divergence makes it numerically difficult to detect the 
vanishing magnetization, in particular for large values of 
the spin. We presume that the magnetization-driven sce- 
nario leads to a weak first order transition to a possible 
intermediate phase and that the velocity- driven scenario 
indicates a strong first order transition to the Neel phase. 

Second, we discussed the importance of a biquadratic 
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0.5 1 1.5 0.5 1 1.5 




ig)E = 144 ± 15 meV, fcc = 5 ih)E = 175 ± 15 mcV, fcc = 5.2 
for 5 = 2. The parameters are given in Tab. IIIII 

exchange. The motivation is the very distinctive spatial 
anisotropy of the dispersion in the pnictides which can- 
not be induced by the small distortive anisotropy. In or- 
der to be able to treat biquadratic exchange on the same 
footing as we treated nearest-neighbor Heisenberg Hamil- 
tonians before we studied two mean-field approaches to 
biquadratic exchange. One is based on the Dyson-Maleev 
representation, the other on the Schwinger boson repre- 
sentation. Surprisingly, even at zero temperature these 
mean-field theories provide distinctively different results. 
We tested both approaches for the Neel phase of the 
nearest-neighbor Heisenberg on the square lattice against 
series expansion data and exact diagonalization. We es- 
tablished that the Dyson-Maleev mean-field theory pro- 
vides reliable results within about 10% of the contribu- 
tion of the biquadratic exchange. 

Third, we applied the developed mean-field approach 
to the columnar phase of the Ji-Jhq-J2-Jc Heisenberg 
model in three dimensions. The aim was to explain the 
magnetic dispersions observed in the undoped iron pnic- 
tides, parent compounds of a novel class of superconduc- 
tors. It turned out that a biquadratic exchange indeed 
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enhances the spatial anisotropy of the spin-wave veloc- 
ities. But the effect is reduced by quantum corrections 
so that it is not large enough for S* = 1 to match the 
experimental situation. Only for larger values of S an 
almost vanishing effective coupling in b direction along 
the stripes of parallel spins is possible. For S = 3/2, 
however, the required exchange couplings are unlikely. 
For instance, the biquadratic exchange would have to 
be larger than the nearest neighbor exchange. But for 
S = 2 perfect agreement of the dispersion is obtained 
without invoking any spatial anisotropy of the model it- 
self. Note that S = 2 corresponds to the chemical valence 
Fe^"*" which follows from an ionic balance of charge in the 
pnictides. This valence implies four holes in the d-shell 
so that the maximum spin according to Hund's rule is 
S = 2. We provided also results for the dynamic struc- 
ture factor for the matching model at = 2 and for 
S — 1 for comparison. 

Of course, there are open points which are beyond the 
scope of the present article. First, there is the relatively 
small staggered moments seen in experiment which indi- 
cate rather smaller than larger spin values. The compli- 
cated local electronic situation and the residual itineracy 
of the charges are likely candidates to cause this discrep- 
ancy. But we mention that recent experiments also found 
large moments in related substances which are also su- 
perconducting. Second, the spin-only model considered 
here does not include the important charge degrees of 
freedom which generate Landau damping and the bad 
metallic behavior of the externally undoped parent pnic- 
tides. So life-time effects are not treated. 

Further research is called for: On the level of local- 
ized spin models the investigation of the recently found 



high-spin substances appears to be interesting but may 
require the more challenging treatment of canted mag- 
netic order. On the level of itinerant charges the influ- 
ence of doped charges has to be studied. These charges 
can be externally doped or internally self-doped between 
the bands. Furthermore, it would be interesting to learn 
more about the values for biquadratic exchanges from 
density- functional calculations. 
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Appendix A: Dyson-Maleev representation 

1. Antiferromagnetic exchange 

For antiparallel orientation of coupled spins, the 
Hamiltonian ([5]) is given in Dyson-Maleev representation 
by 

with 



rj afm 



Tjafm 



jY,i-S' + s[n,+ n, + bib] + - - 1 {bln^b] + k,n,b^) | 



(A2a) 



, fii -f n. + Iq 



+ ^n'f +n^j+ ^nSj + Iq + i.fi, + fij) Iq + Iq {n, + fij) + 



S 



2 (hifi^ + n^fij) + nifijla + lofiihj + - {1^1^ + I^Ig + [hi + fij) /„ + /g {fii + hj)) 



(A2b) 



where 



h -.^bhl + b.b 



In ■■=b]h,bl+b,n,,b^ 



(A3a) 
(A3b) 



■ "I ■"I "J ' "I 

with i,j being the two coupled sites. For the mean- field 
decoupling the reader is referred to Eq. 



2. Ferromagnetic exchange 

For parallel orientation of coupled spins, the Hamilto- 
nian ^ in Dyson-Maleev representation reads 



rrfm Trim , rrfm 



(A4) 
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where 



(A5a) 



rrfm _ 



(A5b) 



r 



with 



btb 



■I" J I "J "I 



Fo 



b% 



bln^b,. 



(A6a) 
(A6b) 



For the mean-field decoupling the reader is referred to 
Eq. (1251). 
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